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Simulations of the structure and dynamics of fluid films confined to a thick- 
ness of a few molecular diameters are described. Confining walls introduce 
layering and in-plane order in the adjacent fluid. The latter is essential to 
transfer of shear stress. As the film thickness is decreased, by increasing pres- 

X 

$— i ' sure or decreasing the number of molecular layers, the entire film may undergo 

a phase transition. Spherical molecules tend to crystallize, while short chain 
molecules enter a glassy state with strong local orientational and translational 
order. These phase transitions lead to dramatic changes in the response of 
the film to imposed shear velocities v. Spherical molecules show an abrupt 
transition from Newtonian response to a yield stress as they crystallize. Chain 
molecules exhibit a continuously growing regime of non-Newtonian behavior 
where the shear viscosity drops as v ~ 2 / 3 at constant normal load. The same 
power law is found for a wide range of parameters, and extends to lower and 
lower velocities as a glass transition is approached. Once in the glassy state, 



chain molecules exhibit a finite yield stress. Shear may occur either within 
the film or at the film/wall interface. Interfacial shear dominates when films 
become glassy and when the film viscosity is increased by increasing the chain 
length. 

PACS numbers: 68.15.+e, 81.40.Pq, 64.70.Pf, 62.10,+s 
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I. INTRODUCTION 



The rheological behavior of fluids in highly confined geometries is a subject of immense 
technological importance. Knowledge of the behavior of confined fluids is crucial to un- 
derstanding flows and phase transitions in porous media, the stability and dynamics of 
coatings, and friction and wear in boundary lubrication and nanomachines. Our knowledge 
of such thin films has increased greatly in recent years due to the development of faster 
computers and new experimental techniques. Examples of the latter include atomic force 
microscopy the quartz crystal microbalance ||, and the surface forces apparatus (SFA) 

In this paper we use molecular dynamics simulations to study the structure and dynamics 
of fluids in a geometry that models the SFA P,|5|,|6|J^,P,p|,p!0| . In the SFA, fluids are confined 
between mica plates that are atomically flat over contact diameters of order 100/im. The 
separation h between the plates can be varied with angstrom resolution from contact to 
many nanometers, and is monitored using optical interferometry. Normal and shear forces 
on the film are measured as the plate separation is decreased, or as the plates are sheared 
past each other at velocity v. 

SFA experiments have focused primarily on films of simple hydrocarbon liquids, such as 
cyclohexane, dodecane, tetradecane, hexadecane, and the silicone-liquid octamethlycyclote- 
trasiloxane (OMCTS) Cyclohexane and OMCTS are roughly spher- 



ical molecules, while dodecane, tetradecane, and hexadecane are short, symmetric chains. 
Other fluids that have been studied include the branched molecules squalane and isoparaf- 
fin 2-methyloctadecane, and polymers, such as polydimethylsiloxane (PDMS) and perflu- 
oropolyether (PFPE) fT^filjj . The latter is commonly used in the computer industry as 
a lubricant for thin film magnetic disks where operating conditions typically require film 
thicknesses of order 1 — 5nm. 

These experiments have shown that the static and dynamic properties of fluids in highly 
confined geometries can be remarkably different from bulk properties. In most cases, bulk 



behavior persists until separations of order ten molecular diameters. At smaller separations 
there are pronounced oscillations in the normal load as a function of film thickness fli"5| , |16| . 



These are evidence of strong layering within the film [|IT| , |T8|| . Changes in the dynamic 
response are even more dramatic. The effective viscosities fi of thin films rise more than five 
orders of magnitude above bulk values. Relaxation times are fractions of a second rather 
than nanoseconds. When films are sheared rapidly, a wide variety of molecules exhibit the 
same power law decrease of viscosity with velocity, [i oc v~ 2 ^ 3 P, |iO| , [T2| . In some cases, the 
response of thin films becomes solid-like: shear stresses do not relax to zero and there is a 
substantial yield stress [||J|[7||§. Even those long-chain and branched molecules which do 
not exhibit layering, still develop a yield stress after long equilibration or shearing. When 
the yield stress is exceeded, "solid" films often exhibit oscillatory stick-slip dynamics 

Computer simulations have proved to be an effective tool in understanding the origin of 
some of these experimental observations. In addition to revealing the layering that was in- 
ferred from experiment and analytic work |T^,|T^,|1^,|T^,^,^1|] , simulations show that the peri- 
odic potential of the wall lattice induces in-plane order in the adjacent fluid P2|,| 



OOHt 



The degree of in-plane order was found to correlate with the viscosity at the solid/fluid 
interface and thus to determine the boundary condition for fluid flow [ 23| , P7| | . As the film 
thickness decreases, confinement can induce a phase transition to a crystalline or glassy 



phase P5] , |^ , |^ , |^ , |^ , |3^ , |33| , |34| . Near the glass transition, films exhibit a non- Newtonian 



response with the same power law viscosity seen in experiments P, p!0| , |^ , |^ , |31] , |32| , |35[| . Stud- 
ies of crystalline films of spherical molecules showed that stick-slip motion resulted from 



periodic phase transitions between sliding fluid and static solid phases p8|j3T[] . This sug- 
gested that the pervasive phenomenon of stick-slip motion is generally associated with phase 
transitions between distinct sliding and static states. 

In this paper we take a deeper look at the effect of molecular geometry on the structure 
and dynamics of thin films. We present results for both spherical molecules and freely jointed, 
linear chain molecules of varying length n. We first describe the changes in equilibrium 
structure and diffusion that occur as films enter a glassy or crystalline state. The transition 



pressure and temperature are shifted from bulk values by confinement. In some cases the 
nature of the transition is also different than in the bulk, i.e. the system enters a glassy 
rather than crystalline state. We next consider the dynamic response to shear. Changing 
the wall/fluid interaction parameters leads to very different types of flow profile: shear at 
the wall/fluid interface, uniform shear within the film, and shear between layers that are 
strongly adsorbed onto each wall. However, in each case we find the same power law drop 
in viscosity with velocity at constant normal load, \x oc v~ 2 / 3 . Simulations at constant film 
thickness produce a less rapid fall in viscosity. Data for the diffusion and relaxation times 
are fit to the free volume theory for glass transitions [06 ]. Results for slower velocities and 
higher pressures than studied previously p9|j30|j3l|j3^[ , reveal that the glassy phase exhibits 
a yield stress (implying \x oc v^ 1 ). This has also been observed in recent experiments [[37j . 



The organization of the rest of the paper is as follows. The next section describes the 
simulation techniques used in our study. We then present results that illustrate the effect 
of confinement on the equilibrium properties of spherical and chain molecules. Section IV 
examines the corresponding changes in the dynamic response of films to an imposed shear 
velocity. Section V contains a brief summary and discussion. 



II. SIMULATION GEOMETRY AND METHOD 



In order to mimic the SFA, our simulations were performed in a planar, Couette geometry 
(Fig. [l]). A thin film of spherical or short-chain molecules was confined between two solid 
plates parallel to the xy plane. Edge effects were minimized by applying periodic boundary 
conditions within this plane. Each wall consisted of 2N W atoms forming two [111] planes 
of an fee lattice with nearest neighbor spacing d. Since mica is much more rigid than the 
confined films, we simplified the simulations described below by fixing wall atoms to their 
lattice sites. Previous work |23] shows that allowing the wall atoms to move does not produce 
qualitative changes in the behavior of the film. The main effect is to produce a small decrease 
in the ability of the wall to order the adjacent fluid at a given temperature T and pressure 
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P, and a corresponding decrease in interfacial viscosity. The definition of the film thickness 
h becomes ambiguous at the atomic level. As shown in the figure and discussed below, h 
was defined to exclude the volume occupied by wall atoms. 

The chain molecules comprising the fluid were simulated using a well-studied bead-spring 
model []2"§U5g|] . Monomers separated by distance r interacted through a truncated Lennard- 
Jones (LJ) potential, 

4ef(o-/r) 12 - (a/r) 6 l , r < r c 
V LJ (r) = 1 LV ' ' V ' ' 1 (2.1) 

[ r > r c , 

characterized by energy and length scales e and a, respectively. In some simulations the in- 
teraction was made purely repulsive by setting the cutoff radius r c /o = 2 1 / 6 . This minimized 
the computation time and the temperature dependence. In other simulations r c /a = 2.2. 
Increasing r c decreases the pressure needed to reach a given density, and increases the vis- 
cosity at that density. These changes do not alter the nature of the transitions in structure 
and dynamics that we describe below. However, they do shift the location of transitions in 
the (P, T) phase diagram. 

The number n of monomers in a chain was varied from 1 (spherical molecules) to 24. 
Adjacent monomers along each chain were coupled through an additional strongly attractive, 
anharmonic potential, 

vCH{r) \-\kRlHl-{r/R ?],r<R (22) 
oo , r > R Q . 

with R = 1.5a and k = 30ea~ 2 . These values of R Q and k have been used in previous 
studies p9|j38[| of polymer melts at comparable density and temperature. They were shown 



to eliminate unphysical bond crossings and bond-breaking. Another important aspect of 
this choice of parameters was that k was not so strong that there was a separation of time 
scales between motions induced by V LJ and V . 

Wall and fluid atoms also interacted with a LJ potential, but with different energy and 
length scales: e w , a w and r cw . Increasing the ratio e w /e, increases the viscosity of the 



wall/fluid interface relative to that within the fluid. The effective viscosity at the interface 



is also strongly affected by the relative size of wall and fluid atoms ||23|| . When the sizes 
are equal, it is easy for fluid atoms to lock into epitaxial registry with the substrate, and 
the viscous coupling is maximized. When the sizes are mismatched, the coupling is weaker. 
Chain molecules have two characteristic sizes. The intermolecular potential V LJ has a 
minimum at 2 1 / 6 cr, while V CH reduces the intramolecular separation to about 0.96a in our 
simulations. The presence of an extra length scale frustrates epitaxial order as discussed 
below. 

The simulations were performed in an ensemble where the number of monomers Nf, the 
system temperature T = l.le/ks, an d the normal pressure or load P± exerted on the top 
wall were all constant. The latter constraint allowed the plate separation h to vary, as in 
experiments. To maintain constant load in the z direction, we added the following equation 
of motion for the top wall: 

Z=(P ZZ -P ± )A/M, (2.3) 

where Z is the average z coordinate of atoms in the top wall, P zz is the zz-component of 
the instantaneous, microscopic pressure-stress tensor at the wall, and A is the area of 
the wall. The mass M is an adjustable parameter that controls the oscillation frequency 
of the wall. A small M leads to unphysically rapid oscillations, while a large M leads to a 
slow exploration of the parameter space associated with h. The times associated with the 
actual mass of mica sheets in an SFA are prohibitively long. We found that our results were 
relatively insensitive to M as long as the period of height oscillations was longer than that 
of the longest phonon period in our finite fluid films. This condition was satisfied in the 
simulations described below by setting M equal to twice the product of the monomer mass 
and the number of wall atoms. 

Shear was imposed by pulling the top plate at a fixed velocity v. The shear stress was 
obtained from two independent and consistent methods: direct calculation of the forces 
exerted on the walls by fluid molecules, or from the xz component of the the microscopic 



pressure-stress tensor [|39| . In the SFA, the transverse alignment of walls in the y direction is 
not rigidly fixed. The top wall can displace in this direction as it slides in response to stress 
from the film. We incorporated this motion in most of the simulations described below by 
coupling the y degree of freedom of the top plate to a spring of force constant k v . The extra 
equation of motion is: 

Y =(P yz A-K y Y)/M, (2.4) 

where Y is the average y coordinate of atoms in the top wall. The transverse stress, P yz 
can generally be minimized by displacements of order d. We used a value of K y /M = 
Omem^cr- 2 that was small enough to allow such displacements, and large enough to pro- 
hibit displacements that were comparable to the system size. Simulations with variable Y 
gave slightly smaller shear stresses than those with fixed Y, but did not change the scaling 
of shear stress with shear rate. 

Constant temperature was maintained by coupling the y component of the velocity to a 



thermal reservoir [i0|. Langevin noise and frictional terms were added to the equation of 
motion in the y direction, 

my = -^-(V LJ + V CH ) -mTy + W, (2.5) 
dy 

where T is the friction constant that controls the rate of heat exchange with the reservoir, 
and W is the Gaussian distributed random force from the heat bath acting on each monomer. 
Note that the variance of W is related to T via the fluctuation-dissipation theorem, 

< W{t)W{t r ) >= 2mk B TT5(t - if) . (2.6) 

If T is too small, the energy dissipated in shearing the system will cause the temperature to 
rise. If T is too large, the random force is not integrated accurately. We found the system 
was well-behaved with V = 2r~ 1 , where r = (ma 2 / e) l l 2 is the characteristic LJ time. The 
equations of motion were integrated using a fifth-order Gear predictor-corrector algorithm 
with a time step At = .005r. 
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The major difference between our theoretical ensemble and experiment is the constraint 
of constant Nf. In experiments, fluid can drain from the compressed region between the mica 
plates to equilibrate with the external chemical potential. However, the equilibration time 
may be extremely long due to the small ratio between h (~ lnm) and the contact diameter 
(~ lOOyum) and the immense viscosities of confined films [|41||. Equilibration is more likely in 



steady sliding experiments where the contact area changes completely flplJlTi , than in the 



experiments with small (~ lOnm) sinusoidal motions that we compare to here □.[Tnj.[TT| 



To determine Nf we used the fact that the film orders into well defined layers of fairly 
constant density p^. We set Nf = miApi where mi was the desired number of layers. This 
guaranteed that the nominal density of films was independent of the number of layers. Equi- 
librium configurations were obtained through a sequence of runs. The load was decreased to 
a small value for 10 4 At, so that chains could move rapidly. The load was then increased se- 
quentially, allowing at least 10 5 At for equilibration. To check that equilibrium was reached, 
we compared runs with different equilibration times and starting configurations. We also 
compared results with increasing and decreasing load, and runs equilibrated at fixed load 
and high temperature. The film thickness h was defined in terms of the distance ho between 
the innermost wall layers as h — homi/(mi + 1). This is equivalent to subtracting a distance 
of half the layer spacing from ho for each wall to account for the volume taken up by wall 
atoms. Other methods of correcting for the size of wall atoms give essentially the same 
results. At most the viscosity is changed by a constant multiplicative factor of order unity. 



III. EQUILIBRIUM STRUCTURE 

Two types of structure are induced in the fluid adjacent to a flat solid surface: Layering 
normal to the surface and epitaxial order in the plane of the surface. When a fluid film is 
confined between two solid surfaces, both types of order may span the entire film. Experi- 



mental consequences are oscillations in the normal force with film thickness [I5|, and phase 
transitions to solid states that resist shear forces j|,|6],|7||| . 
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Layering has been studied extensively for a variety of molecular geometries 
g7J00^g^3^3^1^30g3ilS- lt is induced b y the monomer pair correlation func- 



tion g(r), and the sharp cutoff in fluid density at the wall. Fig. ^ shows plots of the density 
as a function of the distance between walls for monomers and 6-mers. Note the well-defined 
density peaks near the walls, and the decay of oscillations into the center of the film. 

In general, the sharpness and height of the first density peak are determined mainly 
by the wall/fluid interaction, and increase with e w or P±. The rate at which the density 
oscillations decay is determined by the decay of correlations in g(r). Except near critical 
points, the decay is a few molecular diameters. For fluids of simple spherical molecules, 
pronounced density oscillations can extend up to ~ 5a from an isolated solid surface. The 
competition between intra- and inter- molecular spacings leads to less pronounced layering 
with chain molecules (Fig. 0). The degree of layering is fairly independent of chain length 
once n exceeds 6. Simulations with realistic potentials for alkanes show more pronounced 
layering near the wall due to a transition to an extended state of the chains in the first layer 



In-plane order is induced by the in-plane variation, or corrugation, in the potential 
from wall atoms. As for layering, the rate of decay is determined by g(r). It has been 



studied extensively for spherical molecules |[17|Jlq , |^j23| , |25|j26| , but only recently for chain 
molecules |]35| , |46|j47|| . In part, this is because most studies of chain molecules near surfaces 
have neglected the discrete lattice structure of the solid surface. In-plane order plays an 
essential role in determining the dynamics of fluids near solid interfaces. For example, we 
have shown that flow boundary conditions near solids are well correlated with the amount 
of in-plane order, and not with the degree of layering |23|;[27 . 



The amount of in-plane order depends upon the relative size of wall and fluid atoms, as 
well as the strength of their interaction [p^fl . If the spacing between fluid atoms is equal to 
that of wall atoms, the atoms can lock in to all the local minima in the wall potential. In an 
earlier study of simple spherical molecules |23| we found crystallization of one or two fluid 
layers adjacent to a solid interface. If the wall and fluid atoms have different sizes, epitaxial 
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order is frustrated [^j. As we now show, the presence of two different length scales in our 
chain molecules also frustrates in-plane order. 

Two measures were used to quantify the degree of in-plane ordering. The first was the 
spatial probability distribution pi(x, y) of monomers in the Zth. layer relative to the unit cell 
of the solid lattice. The second was the two-dimensional static structure factor S(k x , k y ) 
evaluated in the Ith layer according to 



Si(k) = l/N t 



e 
i=i 



(3.1) 



where iVj was the number of monomers % within the layer. 

Figures ^ and ^ show p{f) and S(k) for a two layer film at two loads. Values of S 
and p for each layer were averaged for 250r. Then results for layers related by symmetry 
about the middle of the film were averaged to improve statistics. At the lowest applied load, 
P± = 4e<r~ 3 , there is little evidence of any in-plane order in p for either n = 1 or 6. The 
structure factor S shows weak rings that are characteristic of liquid order, and small peaks at 
the reciprocal lattice vectors Gi of the solid surface. These peaks represent a linear response 
to the surface potential. Their strength is best expressed P3[ in terms of the Debye- Waller 



factor e~ 2W = S(Gq)/Ni at the smallest reciprocal lattice vector Gq. The Debye- Waller 
factor is unity in an ideal crystal at zero temperature. Thermal fluctuations decrease e~ 2W 
to about 0.6 at the melting point of a bulk crystal | 48| . The value of e~ 2W at P± = 4ea~~ 3 
is only 0.06 for both n = 1 and 6. Such small values are clear evidence of liquid structure. 

As the load increases, the walls induce more in-plane order. At P± = 16ea~ 3 , there are 
well-defined modulations in p(r). The peaks in p are located above gaps in the adjacent solid 
layer. For n = 1, the monomers are nearly always over these gaps. Less order is induced in 
films of chain molecules because of the different intra- and inter-molecular spacings. Note 
the pronounced ridges that connect the peaks in p(r). These follow the lowest energy path 
between the gaps in wall atoms. 

The monomer results for S(k) and p(r) clearly indicate crystalline order at P± = 16eo~~ 3 . 
For example, the mean-squared displacement of a fluid atom from the peaks in p(f) is 
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~ 0.017<i 2 , which is below the Lindemann criterion for melting [^,^]. Furthermore, S(k) 
shows many sharp Bragg peaks, including several higher-order harmonics. The Debye- 
Waller factor, e~ 2W = 0.78, is above the bulk melting criterion and is consistent with 
the mean-squared displacement. In contrast, the chain molecules never crystallize within 
our simulation times. The degree of order increases with P± and then saturates. The 
Debye- Waller factor is only 0.43 in Fig. f|. 

Figure |^ shows the load dependence of film thickness, Debye- Waller factor and in-plane 
diffusion constant D in 2 layer films p9| . Note the sharp transitions in the monomer results 
at P± = 7ecr~ 3 . At lower loads the structure is liquid-like and diffusion is rapid. At higher 
loads the film has crystallized and diffusion is suppressed. The lowest value of e~ 2W within 



the crystalline phase (0.62) is very close to that in bulk crystals at their melting point [48] 
There is a similar phase transition in bulk films, but it occurs at P± = 12e<j~ 3 . Using typical 
values for e = 2007^ and a = 3A, this corresponds to a shift of 0.5 GPa. Note that the 
absolute loads are sensitive to the potential cutoffs r c and r cw , but the shift in the transition 
point is less so. 



Results for chain molecules show more gradual changes |29||. The Debye- Waller factor 



saturates at about 0.4, and remains nearly unchanged if the load is increased by an addi- 
tional order of magnitude. The only evidence of a transition to a new phase comes from 
measurements of dynamic quantities, such as the diffusion constant. The value of D drops 
rapidly below that in a bulk fluid at the same pressure. Studies of the viscous response, 
described in the next section, also show a dramatic slowing of molecular motions. These 
findings are evidence that the film undergoes a glass transition. Extrapolation to D = 
using free volume theory indicates a transition near 

Pf = iQ e a^ (Fig. [10]). At loads above 
this value we find that films exhibit solid behavior over the time scales of our simulations: 
there is no relaxation of applied shear forces and diffusion is undetectably small. Note that 
diffusion in bulk films extrapolates to zero at a much higher pressure, P^ lk ~ 24ecr~ 3 . As 
for spherical molecules, confinement produces large shifts (~ 0.8GPa) in the bulk transition 
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pressure. These shifts are comparable to the glass transition pressures (~0.5GPa) of bulk 
lubricants at room temperature p0]| , and one may expect that many lubricants vitrify in 
thin films. 

Although chain molecules do not crystallize, there are pronounced changes in their con- 
figuration and orientation as they approach the glass transition. Figure |^ shows the dis- 
tribution of end-to-end distances R\ n for n = 6 at three different pressures. At the lowest 
load, P± = leer -3 , the film thickness h = 2.94a is larger than the radius of gyration of a 
free three-dimensional chain [ 5l| , and the density is nearly independent of z. As a result the 



distribution of end-to-end distances is close to that for an ideal Gaussian chain [p2j| . The 
major difference is a cutoff in the distribution at small separations due to the hard core 
repulsion between monomers. As P±_ increases, the distribution becomes more structured, 
indicating that monomers are being locked into a discrete set of locations. One also sees that 
highly stretched configurations become less likely. Previous studies of polymers confined to 
two dimensions show that they maximize their entropy when each polymer segregates into 
its own region of the plane EBJ . The result is a collapse into a compact configuration of the 



chain. As P±_ increases from 1 to 16ecr 3 , h drops from a value greater than the radius of 
gyration to a smaller value, and the behavior crosses over from three- to two-dimensional 

Figure [7| shows how the orientation of intramolecular bonds varies with load. We use the 
conventional coordinate system where 9 is the polar angle relative to the z axis and is the 
azimuthal angle within the xy plane and relative to the x axis. Since the sequence along 
the chain is arbitrary, the sign of the vector connecting nearest neighbors is not defined. 
Thus we calculated the distribution of cos 2 9 for all intra-molecular bonds. The distribution 
for completely random orientations, 1/(2 cos 9), is shown by dashed lines in the plots. At 

= 4ecr~ 3 , the distribution is nearly random. For P± = 16eo~~ 3 one finds sharp peaks at 
preferred orientations. The suppression of other orientations limits the ability of molecules 
to realign and results in a decreased diffusion rate. At the preferred values of cos 2 9 = and 
0.6, bonds connect monomers in the same (9 = 90°) or adjacent (9 ~ 40°) layers. Similar 
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peaks appear in the distribution of cos 2 <fi near and 0.75, or = 90 and 30°. These are the 
directions of bonds between monomers that have locked epitaxially into minima in the wall 
potential or lie along the density ridges in Fig. [|. 

The structural changes in 2 layer films of chains with n > 6 are nearly the same as the 
6-mer results just described. In general, the structure becomes insensitive to chain length 
when the film thickness is much smaller than the radius of gyration in a bulk melt. As we 
now show, the shear response of the system may also become insensitive to n in this limit. 



IV. RESPONSE TO STEADY SHEAR 



SEA experiments probe the dynamic response of films by measuring the shear force 



on the walls as a function of velocity P,ffO,B3. This measurement necessarily combines 



information about the shear response of the film and the interface. There is no way to 
determine whether shear occurs primarily within the film or is localized at the interface. It 
has also been impossible to examine structural changes in shearing films. 

Experimentalists generally assume a no-slip boundary condition at the wall/fluid inter- 
face, i.e. that all shear occurs within the film. The shear rate 7 = dv x /dz is then given by 
v/h. An effective viscosity of the film fi can be determined using the macroscopic relation 
between force and shear rate in a lubricant film 

/ = A*7 = t*v/h (4.1) 

where / is the force per unit area or shear stress. In a simple Newtonian fluid, /x is indepen- 
dent of shear rate and / rises linearly with v. In most non-Newtonian fluids, \i decreases 
with v, and / rises sublinearly |]55| ,j56||. 

Typical force- velocity curves for chain molecules are shown in Fig. |8], where n = 6, e w = e 
and mi = 2. The top wall was sheared at constant velocity v in the x direction, and allowed 
to adjust its registry in the y direction in response to shear stresses as described above. 
Results for P± = 8 and 12e<7~ 3 rise rapidly from zero and then saturate as v increases. The 
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film is highly non-Newtonian and, as shown below, / scales roughly as v 1//3 . At P± = 16 
and 20ecr~ 3 , the force approaches a constant value or yield stress as v — > 0. This behavior 



is characteristic of solid-on-solid friction |)7|], and is further evidence that the film is in a 
glassy state at these loads. 



Granick and coworkers |J|nj have plotted their shear response curves in terms of the 
effective viscosity and shear rate defined in Eq. |4.1| . Fig. |^(a) shows the data of Fig. |8| 
replotted in this manner. A Newtonian regime with constant viscosity /io is seen at the 
lowest loads and shear rates. As 7 increases, the system can no longer respond rapidly 
enough to keep up with the sliding walls. When P± is below the glass transition load 
there is well-defined crossover shear rate j c at which the viscosity begins to drop. The 
shear-thinning appears to follow a universal power law \i oc v ~ 2 ^ 3 (implying / oc v 1 ' 3 ) that 
is indicated by a dashed line on the figure P, |l0| , |29 . 



The inverse of j c corresponds to the longest structural relaxation time of the film |55| , |56 
As P± rises toward P±, ^ c drops to zero and the relaxation time diverges. This slowing of 
dynamics in the film is directly related to the decrease in diffusion constant with P± seen in 
Fig. |5|, and the rapid rise in /x with P± seen in Fig. |9|. 



One common description of glass transitions is the free volume model It states that 



7 C and D should vanish as exp(—ho/(h — h c )) where h c is the film thickness at the glass 



transition |5T|. Fig. [10| shows fits of both quantities to this form with h c corresponding to 
a glass transition near P± = 16e<r~ 3 . There is considerable debate about the proper scaling 
of dynamics near a glass transition, and even about the existence of a sharp transition. The 



success of the fit in Fig. |10| should be taken as evidence that the transition we see is like 
other glass transitions, rather than evidence for the free volume theory. 

When Pj_ exceeds P^, there is a qualitative change in the response at small 7. As shown 
in Fig. ^ /i falls more steeply with v at P±_ — 16 and 20ecr~ 3 . The slope on this log-log plot 
approaches —1 as 7 decreases. This implies a constant shear force or yield stress, as already 
seen for these loads in Fig. || A constant shear force is typical of solid-on-solid sliding j57| . 
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We show below that the film is nearly rigid at these loads, and that shear is largely confined 



to the wall/film interface. This regime was not evident in previous simulations p9| , |30| , |3T| , |32 
because they did not extend to high enough loads and low enough velocities. 

Granick et al. have found a very similar series of changes with load |],[n],[J7|] . Indeed they 
are even more dramatic, because a larger range of time scales is accessible to experiments. 
Increasing the load produced a decrease in 7 C by at least 10 orders of magnitude, and a rise 
in the Newtonian viscosity /io by more than 5 orders of magnitude. At loads where 7 C was 
still non-zero ||, there was a pronounced non- Newtonian regime where n dropped as v ~ 2 / 3 . 
At higher loads there is a crossover to constant shear force |[37|| . 

A decrease in viscosity with increasing shear-rate is normally accompanied by structural 
changes that reflect the fact that the system no longer has time to relax back to the equi- 
librium state |55|,56| . The main structural change in our constant load simulations is an 



increase in the film thickness with shear rate. Shear produces an extra normal force that 
separates the walls. This facilitates sliding and lowers fi. Note in Fig. |9|(b) that h is constant 
in the Newtonian regime for a given load and begins to rise only when 7 exceeds the value 
of 7 C . Once the film is in the glassy phase, changes in h extend to the lowest practical shear 
rates. The changes in h may be coupled to structural changes within the film, but these are 
difficult to detect until relatively high shear rates (7 > O.lr -1 ). In this regime, the number 



of intramolecular bonds which connect layers begins to decrease to facilitate flow |35l , f46| , |47 

The above results imply that fi should drop less rapidly with 7 if h is fixed. This is 
verified by the constant h data shown in Figure O. We find that fi decreases with a power 



law near -0.5 in this ensemble p9| , |35]| . Previous experiments have been done at fixed P±_ 
and it would be interesting to test our simulation results against experiments at fixed h. 
Studies of structural changes in our fixed h simulations show that the layers of monomers 
move away from the walls to facilitate interfacial shear and the layering at the wall becomes 
slightly more defined. As above, there is little structural change within the films. 



Figure O shows how / depends on chain length in 2 layer films at = lQea 3 . Note 



that the results become independent of chain length for large n. Only the results for n = 1 
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and 3 are noticeably different. This length independence may seem surprising since the 
bulk viscosity rises rapidly with n. However, at these high loads, films of long chains have 
frozen into a glassy state and all the shear occurs at the interface. As noted in the previous 
section, the structure of the interface and the resulting shear coupling are insensitive to n 
when chains are long enough to span the system. 

In Fig. fi and 7 were obtained from / and v by assuming a no-slip flow bound- 
ary condition (Eq. [4.1|) . This assumption is invalid when shear localizes at the interface. 
Mechanical equilibrium requires that the shear stress be uniform throughout the system, 



however the shear rate may vary with position. In Fig. [13] we show the variation of v x with 
z for a number of systems. The actual shear rate within the film, 7/jj m = dv x /dz, is given 
by the slope of the flow profiles. For the monomer results shown in Fig. ]13|(a), there is 
relatively little shear at the walls and A ffu m ~ 7. As n increases to 6 (Fig. [13] (a)), shear 
becomes predominantly confined to the wall/film interface. Fig. [14] shows the variation of 
ifUm/i with n at P± = 16ea~ 3 . This ratio would be unity if a no-slip boundary condition 
applied. While the no-slip approximation is very accurate for monomers, it fails rapidly as 
n increases to 2 and beyond. For long chains, only 15% of the shear occurs within the film 
at P± = 16ecr~ 3 . This type of flow profile is frequently called plug-like flow. 

The power law shear-thinning shown for two layer films in Fig. || is quite universal. Fig. 



15 shows results for a system where e w was increased to 3e in order to decrease the amount of 



shear at the interface [29]. Note that the glass transition can be approached in two different 



ways because the same load can be obtained with different numbers of layers fL5| . Panel 



(a) shows results at fixed load with a decreasing number of layers, and panel (b) shows the 
effect of increasing load at mi = 4. In either case, there is an increase in the Newtonian 
viscosity and a decrease in j c as the fluid becomes more confined. The shear-thinning region 
shows the same power law dependence seen in Fig. || and experiment |],|llj] . If the load is 
increased further into the glassy regime, shear will localize at the interface. In this limit one 
finds a regime where the shear force is velocity independent. 

A typical flow profile for the system just discussed is shown in Fig. 0(b). Note that 
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while all the shear occurs within the film, it is not spread uniformly throughout the film. 
Over the entire range of power law shear thinning, shear is localized at a plane near the 
center of the film. Closer examination of particle motions reveals that each molecule adsorbs 
tightly to one of the two walls. These adsorbed films then slide past each other. Hence, 
the power law shear thinning is still associated with interfacial sliding. The only difference 
from the system of Fig. [5] is that the interface is within the film rather than at the wall film 
interface. 

To determine whether uniform shear could produce the same power law response, we 
examined several other sets of parameters. Fig. |13](b) also shows a flow profile for e w /e = 1, 
P± = 8ecr -3 , and r cw = r c = 2.2a. Note that the shear rate is constant within the film 
(v x is linear). Roughly half of the shear occurs in the film and the other half at the two 
wall/film interfaces. Fig. |16| shows the variation of \i with 7. As for other systems, there is 
a substantial range where [i oc 7~ 2 / 3 . Thus we conclude that even films with uniform shear 
exhibit the same power law near the glass transition. Note that if the interface and film 
had different power laws the flow profile would change with 7. We find the same division of 
shear between interface and film over the entire range of shear rates. 



V. DISCUSSION AND SUMMARY 



The simulations reported here, and other related work, show that confinement can pro- 
duce a rich variety of phenomena in nanometer scale films. Structural changes, including 



layering and in-plane order, are now well-documented []lq, P0| , pi] , P3 ) H1 ■ We have shown here, 
that both types of ordering are smaller in chain molecules due to the presence of more than 
one length scale. This effect should be even more pronounced in films of branched or other 
more complicated molecules. 

Confinement may eventually lead to a phase transition within the film (Fig. [Sp 
2^.BD].pT|,|S2| . In some cases, the transition is just a bulk phase transition that is shifted 



to a lower pressure or higher temperature. This is the case for crystallization of spherical 



molecules by walls whose atoms have nearly the same size. In other cases, the nature of 
the transition is changed. Chain molecules tend to be trapped in a glassy state. Spherical 
molecules can form a glassy phase or a new crystal structure when the walls are amorphous 
or epitaxy is frustrated due to size differences |23| , p8| . Glassy phases also occur when the 
crystalline axes of the confining walls are not aligned, because the walls induce contradictory 



ordering tendencies on the film [p8| . The mica walls in the SFA are generally not aligned 
and glassy phases may be the rule in experimental studies. This would explain why the 
same power-law non-Newtonian behavior is observed in films of chain molecules, like alka- 
nes, and more spherical molecules, like OMCTS ||10]. An interesting observation from our 



simulations is that thicker films (3 or 4 layers) are more able to accommodate the conflicting 
influences of the walls and may have a higher yield stress than films with only 1 or 2 layers. 
Such films typically melt when they begin to slide and will be the topic of a later paper. 

Three distinct types of dynamic shear response were found. Newtonian behavior occurs 
at large thicknesses and temperatures and at low pressures. In the opposite limits there is 
a constant shear stress or yield stress as v — > 0. The transition between the two types of 
behavior occurs abruptly if the film crystallizes, and gradually if it vitrifies. The third type 
of response occurs near the glass transition where there is an extended range of power law 
shear thinning /i oc 7~ 2 / 3 (/ oc v 1 ' 3 ). The power law holds from a crossover shear rate j c 
up to the frequencies of typical phonons. As the glass transition is approached, 7 C drops 
rapidly to zero. There is a corresponding drop in the diffusion constant and a rise in the 
value of /i at 7 < j c . 

These types of response describe the behavior of the entire system, including film and 
wall/film interfaces. Although the shear stress is independent of position, the shear rates 
within the film and at the interface may be very different. The shear rate at the interface 
decreases as the strength of the wall fluid coupling (e w and/or r cw ) increases, and when d is 
such that monomers can easily lock into epitaxy with the wall. The shear rate within the 
film depends strongly on n, mi and their relative sizes. When chains are strongly adsorbed, 
shear may become localized at the midpoint of the film. 
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All flow profiles seem to be associated with the same power law shear thinning. A 
different power law was obtained only when the ensemble was changed from constant load 
to constant thickness. For constant film thickness jx decreases as 7 _1//2 . Recent simulation 
studies of the bulk shear viscosity for the same molecular potential find a very similar shear 
thinning exponent at constant volume . It would be interesting to compare these results 
to bulk simulations at fixed pressure. 

Several models have been developed to explain the power law shear thinning observed in 
experiment and simulations |60| , |6llj62]| . All the models find shear thinning with an exponent 
of 2/3 in certain limits. However, they start from very different sets of assumptions and it 
is not clear if any of these correspond to our systems. Two of the models work at constant 



film thickness |60|j6T| where we find an exponent of 1/2. One of these |61j also finds that 
the exponent depends on the velocity profile, while our results do not. The final model 
|62|| is based on scaling results for the stretching of polymers under shear. While it may be 



appropriate for thick films, it can not describe the behavior of films which exhibit plug-like 
flow. It remains to be seen if the 2/3 exponent has a single explanation or arises from 
different mechanisms in different limits. 

Shear is known to stretch and align bulk polymers. One might expect similar changes 
in thin films, but confinement limits the ability of polymers to rotate in the xz plane. It 
has been suggested that chains might prefer to align normal to the xz plane in this case 
||. We found no tendency for shear-induced alignment in any direction when the film was 
much thinner than the bulk radius of gyration and the shear rate was less than O.lr -1 . Fig. 



17| shows the decay of order in a film where all molecules were originally aligned in the x 
direction. The projections of the vector between the ends of the chains on to the x and y 
axes are plotted against time. There is a rapid decrease in anisotropy over ~ lOr, and the 
ratio of the projections (Fig. |rT|(c)) reaches unity after about 450r. There does appear to be 
memory of the direction of sliding for even longer times, but this must be reflected in more 
subtle details of the structure, such as the relative positions of monomers in the xy plane. 
Studies of the development of orientational order as the degree of confinement is reduced 
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are underway |k| . 
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FIGURES 

FIG. 1. Projection of particle positions onto the x—z plane illustrating the simulation geometry. 
Periodic boundary conditions are imposed in the x and y directions. Gray circles connected by 
lines show the short chain molecules that make up the film of thickness h. Solid circles indicate 
wall atoms. They sit on the sites of fee [111] layers in the x — y plane. A constant load P± is 
applied in the direction normal to the top wall, and the bottom wall is held fixed. Shear may be 
applied by moving the top wall at fixed velocity v in the x direction. In most simulations, the wall 
was allowed to optimize its registry in the y direction as sliding occurred (Eq. ^4| ) 

FIG. 2. Profiles of the monomer number density p as a function of z at v = for films with 
(a) n = 1, h = 3.41cr and (b) n = 6, h = 3.22cr. Other parameters were e w /e = 1, r cw = r c = 2.2a, 
mi = 4, P± = 8ea~ 3 , ksT = Lie, Nj = 288, N w = 72, and d = 1.2a. The number of monomers in 
bins of width h/80 along the z direction was averaged over 150r. Note that z is normalized by the 
average plate separation, h. 

FIG. 3. Values of p(f) (right) and S(k) (left) for a two layer monomer film at Pj_ = 4 (bottom) 
and 16eer~ 3 (top). The scale for each quantity is adjusted to fill the figure. All data were averaged 
over 250r with v = 0, e w /e = 1, r cw = r c = 2 1 / 6 cr, N w = 144, N f = 288, k B T = Lie, and d = 1.2a. 
Wall atoms lie on the sites of a triangular lattice and are at the corners and center of the unit 
cell shown for p. Sharp Bragg peaks are found in S(k) at the reciprocal lattice vectors of the wall 
surface. At P± = 4ecr~ 3 , the peaks are barely above the circular ridge that is characteristic of fluid 
structure. 

FIG. 4. Results for 6-mers with the same parameters as in ||. Note the pronounced ridges 
between peaks in p for P± = 16ecr~ 3 . 
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FIG. 5. Effect of varying load on (a) plate separation, (b) Debye- Waller factor, and (c) diffusion 
constant for motion parallel to the walls. All data is for v = with e w /e = 1, r cw = r c = 2 1 / 6 <7, 
d = 1.2o", ksT = Lie, N w = 144, Nf = 288, and mi = 2. The open and filled circles denote films 
with n = 1 and 6, respectively. Stars indicate the bulk diffusion constant for n = 6. Results for 
n = 20 (squares) overlap n = 6 results in (a). 

FIG. 6. Effect of varying P± on the distribution T> of end-to-end distances, R\ n , of n = 6 chains 
in systems with v = 0. All parameters as in Fig. [| 

FIG. 7. Effect of P± on chain orientation. Here 6 is the angle between the z-axis and a bond 
between adjacent monomers on a chain. Results for the distribution function T> were averaged over 
250t for the system of Fig. [|. 

FIG. 8. Force per unit wall area, /, as a function of v at the indicated loads. Other parameters 
as in Fig. ||. The force saturates at large velocities. In the limit v — * 0, / goes to when P ± < Pf 
and to a constant when P± > P^. 

FIG. 9. Data of Fig. || replotted in (a) as log/i vs. log 7. The corresponding film thickness at 
each load is shown in (b). At low Pi and 7 the value of fj, is a constant. Above j c , [x begins to 
drop as 7~ 2//3 , and h starts to rise. For P± > P±, ji drops as 7 -1 at low 7. The dashed and dotted 
lines indicate slopes of —2/3 and —1, respectively. 

FIG. 10. Fit of variation in dynamic quantities % and D to free volume theory. The logarithm 
of the inverse of both quantities should scale like ho/(h c — h) where h c is the thickness at P±. In 
the fit, ho = 1.58(7 and h c = 1.57a, implying P^ = 16ecr -3 . 

FIG. 11. Data for the system of Fig. [8|, but with fixed h instead of fixed load. Note that the 
data now indicate that fi falls roughly as j^ 1 / 2 (dashed line) above j c . 

FIG. 12. Frictional force per unit area / as a function of v and n at P± = 16ecr~ 3 . Other 
parameters as in Fig. ||. 
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FIG. 13. Flow profiles for different interaction potentials. Panel (a) shows results for the system 
of Fig. H at v = 0.05<7T _1 with n = 1 (open circles) and n = 6 (closed circles). Stars in (b) are for 
n = 6, mi = 4, P± = 6ea~ 3 , e w /e = 3, r c = 2 1 / 6 <r, r cw = 2.2a, and v = 0.02crr _1 (see also Fig. |l5|) . 
Squares are for n = 6, mi = 4, P± = 8ecr~ 3 , e w /e = 1, r c = 2.2a, r cw = 2.2a, and v = 0.08o"t -1 
(see also Fig. ^). Particle velocities in the x direction, v x , were averaged within layers and over 
at least 250r. These averages were then normalized by the wall velocity v and plotted against the 
center of each layer z divided by h. The flow profile for a no-slip boundary condition is indicated 
by the dotted line in (b). 

FIG. 14. Ratio of the actual shear rate within the film to the value inferred from a no-slip 
boundary condition as a function of n. All other parameters as in Fig. |5[ For n > 3 only about 
15% of the shear is within the film. The rest is at the two interfaces. 

FIG. 15. Plots of jjL vs. 7 at (a) fixed P± = 4e<7 -3 and varying mi, and (b) fixed mi = 4 and 
varying P±. Dashed lines have slope -2/3. Both panels show the same behavior as Fig. ||(a) as 
the glass transition is approached. Other parameters were € w /e = 3, r cw = 2.2a, k B T = Lie, 
r c = 2 1 ' 6 <7, n = 6, and d = a. 

FIG. 16. Plot of fi vs. 7 at fixed P± = 8ecr~ 3 for mi = 4, e w /e = 1, r cw = r c = 2.2a, d = 1.2a, 
and ksT = Lie. As before, the dashed line has slope -2/3. 

FIG. 17. Time dependence of mean-squared projections on to the x and y axes of the vectors 
connecting ends of chains. The initial state was fully aligned along the ^-direction. The ratio of 
the projections decays to unity by the end of the figure, showing that the final state is isotropic. 
Longer runs and runs from other starting states show the same final isotropy. 
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